/*
* Cholesky.cs
* 
* Copyright (c) 2006, dnAnalytics. All rights reserved.
*/

using System;
using dnAnalytics.Math;
using dnAnalytics.Resources;
using dnAnalytics.Exceptions;

namespace dnAnalytics.LinearAlgebra.Decomposition
{
    /// <summary>
    /// Computes the Cholesky factorization of a symmetric, positive definite <see cref="Matrix"/>.
    /// </summary>
    /// <remarks>Cholesky factorization decomposes symmetric, positive definite matrices.
    /// This class does not check whether a matrix is symmetric or not. It only uses the upper
    /// triangle of the matrix to compute the decomposition.</remarks>
    public sealed class Cholesky : Algorithm
    {
        #region Fields

        private ICholesky decomp;
        private readonly int order;

        #endregion Fields

        #region Constructor
        /// <summary>
        /// Constructs an Cholesky object for the given matrix.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <b>null</b>.</exception>
        /// <exception cref="NotSquareMatrixException">If <paramref name="matrix"/> is not a square matrix.</exception>
        public Cholesky(Matrix matrix)
        {
            if (matrix == null)
            {
                throw new ArgumentNullException("matrix", Strings.NullParameterException);
            }

            if (matrix.Rows != matrix.Columns)
            {
                throw new NotSquareMatrixException(Strings.MustBeSquare);
            }

            decomp = new DenseCholesky(matrix);

            order = matrix.Rows;
        }

        #endregion Constructor

        #region Public Methods

        ///<summary>Return a value indicating whether the matrix is positive definite.</summary>
        ///<returns><b>true</b> if the matrix is positive definite; otherwise, <b>false</b>.</returns>
        public bool IsPositiveDefinite()
        {
            return decomp.IsPositiveDefinite();
        }

        ///<summary>Returns the Cholesky factored matrix (lower triangular form).</summary>
        ///<returns>The lower triangular Cholesky factored matrix.</returns>
        ///<exception cref="NotPositiveDefiniteException">If the matrix is not positive definite.</exception>
        public Matrix Factor()
        {
            if (!IsPositiveDefinite())
            {
                throw new NotPositiveDefiniteException();
            }

            return decomp.Factor();
        }

        /// <summary>
        /// Copies the Cholesky factored matrix (lower triangular form) into the result <see cref="Matrix"/>.
        /// </summary>
        /// <param name="result">The <see cref="Matrix"/> store factored matrix into.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="result"/> is <b>null</b>.</exception>
        /// <exception cref="NotSquareMatrixException">If <paramref name="result"/> is not a square matrix.</exception>
        /// <exception cref="NotPositiveDefiniteException">If the matrix is not positive definite.</exception>
        public void Factor(Matrix result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if (result.Rows != order || result.Columns != order)
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }

            if (!IsPositiveDefinite())
            {
                throw new NotPositiveDefiniteException();
            }
            decomp.Factor(result);
        }

        ///<summary>Calculates the determinant of the matrix.</summary>
        ///<returns>The determinant of the matrix.</returns>
        ///<exception cref="NotPositiveDefiniteException">If the matrix is not positive definite.</exception>
        public double Determinant()
        {
            if (!IsPositiveDefinite())
            {
                throw new NotPositiveDefiniteException();
            }

            return decomp.Determinant();
        }

        /// <summary>
        /// Solves a system of linear equations, <b>AX = B</b>.
        /// </summary>
        /// <param name="rightHandSide">The right hand side <see cref="Matrix"/>, <b>B</b>.</param>
        /// <returns>The left hand side <see cref="Matrix"/>, <b>X</b>.</returns>
        /// <exception cref="NotPositiveDefiniteException">If the matrix is not positive definite.</exception>
        /// <exception cref="ArgumentNullException">If <paramref name="rightHandSide"/> is <b>null</b>.</exception>
        /// <exception cref="NotConformableException">If the number of rows of the <paramref name="rightHandSide"/> 
        /// matrix differs from the number of rows of the decomposed matrix.</exception>
        public Matrix Solve(Matrix rightHandSide)
        {
            if (rightHandSide == null)
            {
                throw new ArgumentNullException("rightHandSide", Strings.NullParameterException);
            }

            if (rightHandSide.Rows != this.order)
            {
                throw new NotConformableException(Strings.ParameterNotConformable);
            }

            if (!IsPositiveDefinite())
            {
                throw new NotPositiveDefiniteException();
            }

            return decomp.Solve(rightHandSide);
        }

        /// <summary>
        /// Solves a system of linear equations, <b>AX = B</b>.
        /// </summary>
        /// <param name="rightHandSide">The right hand side <see cref="Matrix"/>, <b>B</b>.</param>
        /// <param name="result">The left hand side <see cref="Matrix"/>, <b>X</b>.</param>
        /// <exception cref="NotPositiveDefiniteException">If the matrix is not positive definite.</exception>
        /// <exception cref="ArgumentNullException">If either <paramref name="rightHandSide"/> or <paramref name="result"/> 
        /// is <b>null</b>.</exception>
        /// <exception cref="NotConformableException"><para>If the number of rows of the <paramref name="rightHandSide"/> 
        /// matrix differs from the number of rows of the decomposed matrix.</para>
        /// <para> - or - </para>
        /// <para>If the <paramref name="rightHandSide"/> matrix has different dimensions than
        /// the <paramref name="result"/> matrix.</para>
        /// </exception>
        public void Solve(Matrix rightHandSide, Matrix result)
        {
            if (rightHandSide == null)
            {
                throw new ArgumentNullException("rightHandSide", Strings.NullParameterException);
            }

            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if (rightHandSide.Rows != this.order)
            {
                throw new NotConformableException(Strings.ParameterNotConformable);
            }

            if (result.Rows != this.order || result.Columns != result.Columns)
            {
                throw new NotConformableException(Strings.ParameterNotConformable);
            }

            if (!IsPositiveDefinite())
            {
                throw new NotPositiveDefiniteException();
            }

            decomp.Solve(rightHandSide, result);
        }

        /// <summary>
        /// Solves a system of linear equations, <b>Ax = b</b>.
        /// </summary>
        /// <param name="rightHandSide">The right hand side vector, <b>b</b>.</param>
        /// <returns>The left hand side <see cref="Vector"/>, <b>x</b>.</returns>
        /// <exception cref="NotPositiveDefiniteException">If the matrix is not positive definite.</exception>
        /// <exception cref="ArgumentNullException">If <paramref name="rightHandSide"/> is <b>null</b>.</exception>
        /// <exception cref="NotConformableException">If the size of <paramref name="rightHandSide"/> does not equal
        /// the number of rows of the decomposed matrix.</exception>
        public Vector Solve(Vector rightHandSide)
        {
            if (rightHandSide == null)
            {
                throw new ArgumentNullException("rightHandSide", Strings.NullParameterException);
            }

            if (rightHandSide.Count != this.order)
            {
                throw new NotConformableException(Strings.ParameterNotConformable);
            }

            if (!IsPositiveDefinite())
            {
                throw new NotPositiveDefiniteException();
            }

            return decomp.Solve(rightHandSide);
        }

        /// <summary>
        /// Solves a system of linear equations, <b>Ax = b</b>.
        /// </summary>
        /// <param name="rightHandSide">The right hand side <see cref="Vector"/>, <b>b</b>.</param>
        /// <param name="result">The left hand side <see cref="Vector"/>, <b>x</b>.</param>
        /// <exception cref="NotPositiveDefiniteException">If the matrix is not positive definite.</exception>
        /// <exception cref="ArgumentNullException">If either <paramref name="rightHandSide"/> or <paramref name="result"/>
        /// is <b>null</b>.</exception>
        /// <exception cref="NotConformableException"><para>If the size of <paramref name="rightHandSide"/> does not equal
        /// the number of rows of the decomposed matrix.</para>
        /// <para> - or - </para>
        /// <para>If the size of <paramref name="result"/> does not equal
        /// the number of columns of the decomposed matrix.</para>
        /// </exception>
        public void Solve(Vector rightHandSide, Vector result)
        {
            if (rightHandSide == null)
            {
                throw new ArgumentNullException("rightHandSide", Strings.NullParameterException);
            }

            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if (rightHandSide.Count != this.order)
            {
                throw new NotConformableException(Strings.ParameterNotConformable);
            }

            if (result.Count != this.order)
            {
                throw new NotConformableException(Strings.ParameterNotConformable);
            }

            if (!IsPositiveDefinite())
            {
                throw new NotPositiveDefiniteException();
            }

            decomp.Solve(rightHandSide, result);
        }


        #endregion Public Methods

        #region Protected Members

        /// <summary>
        /// Delegates the actual computation to the correct Cholesky subtype.
        /// </summary>
        protected override void InternalCompute()
        {
            decomp.Compute();
        }
        #endregion Protected Members
    }
}
